rm(list=ls())
set.seed(654321)

##########
# A Program to Simulated Data from 56 True 

# Load functions to report true Probability y=1

trupry<-function(r,X,Z){
  
  if(r==85) {
    trupry1 = (((X -1)^2 / -10) + .9) + ((((Z + .75)^2 / 2.5) -.225) *((((X -1.250000075)^2 / -3.750000375) + 0.516666675) - (((X -1)^2 / -10) + .9)))
    return(trupry1)}
  
  if (r==82) {
    trupry1 = (((X -1)^2 / 10) + .3) + ((((Z - 1.125)^2 / -1.25) + 1.0125) *((((X -1)^2 / 2.5) + .5) - (((X -1)^2 / 10) + .3)))
    return(trupry1)}
  
  if(r==86) {
    trupry1 = (((X + 0)^2 / 2) + 0.4) + ((((Z + .125)^2 / 1.25) - .0125) *((((X + 0)^2 / 5) + 0.1) - (((X + 0)^2 / 2) + 0.4)))
    return(trupry1)}
  
  if(r==87) {
    trupry1 = ((1-Z)*((((X + 0.2499994)^2 / 2.9999976) + 0.27916675))) + ((Z)*(((X + 0)^2 / 5 + 0.1)))
    return(trupry1)}
  
  if(r==88) {
    trupry1 = (((X - 1.25)^2 / -5) + 0.9125) + ((((Z - 1.125)^2 / -1.25) + 1.0125) *((((X - 1.25)^2 / -2.5) + 0.825) - (((X - 1.25)^2 / -5) + 0.9125)))
    return(trupry1)}
  
  if(r==80) {
    trupry1 = ((1-Z)*((((X -1.249785775)^2 / 2.142245073) + 0.170874979))) + ((Z)*(((X -1.246268657)^2 / 3.731343284 + 0.083746269)))
    return(trupry1)}
  
  if(r==100) {
    trupry1 = (((X -1)^2 / -2) + 0.8) + ((((Z + .125)^2 / 1.25) - .0125) *((((X -1.248502994)^2 / -14.97005988) + 0.90412515) - (((X -1)^2 / -2) + 0.8)))
    return(trupry1)}
  
  if(r==101) {
    trupry1 = (((X + 0)^2 / 10) + 0.2) + ((((Z - 1.125)^2 / -1.25) + 1.0125) *((((X + 0.25030012)^2 / 3.00120048) + 0.37912497) - (((X + 0)^2 / 10) + 0.2)))
    return(trupry1)}
  
  if(r==96) {
    trupry1 = (((X + 2.77556E-16)^2 / -5) + 0.9) + ((((Z + .75)^2 / 2.5) - .225) *((((X + 0)^2 / -1.666666667) + 0.8) - (((X + 2.77556E-16)^2 / -5) + 0.9)))
    return(trupry1)}
  
  if(r==97) {
    trupry1 = (((X -1)^2 / 1.428571429) + 0.2) + ((((Z - 1.75)^2 / -2.5) + 1.225) *((((X - 1.25)^2 / 5) + 0.0875) - (((X -1)^2 / 1.428571429) + 0.2)))
    return(trupry1)}
  
  if(r==111) {
    trupry1 = (((X + 0.25)^2 / -2.5) + 0.725) + ((((Z + .75)^2 / 2.5) - .225) *((((X + 5.55112E-16)^2 / -10) + 0.9) - (((X + 0.25)^2 / -2.5) + 0.725)))
    return(trupry1)}
  
  if(r==112) {
    trupry1 = (((X + 0)^2 / -1.428571429) + 0.8) + ((((Z - 1.75)^2 / -2.5) + 1.225) *((((X + 0.257575758)^2 / -7.575757576) + 0.908757576) - (((X + 0)^2 / -1.428571429) + 0.8)))
    return(trupry1)}
  
  if(r==89) {
    trupry1 = (((X + 0)^2 / 10) + 0.7) + ((((Z + .125)^2 / 1.25) - .0125)*((((X + 0.249625187)^2 / 3.748125937) + 0.083374963) - (((X + 0)^2 / 10) + 0.7)))
    return(trupry1)}
  
  if(r==83) {
    trupry1 = (((X + 5.55112E-16)^2 / -10) + 0.5) + ((((Z - 1.75)^2 / -2.5) + 1.225)*((((X + 0)^2 / -2.5) + 0.9) - (((X + 5.55112E-16)^2 / -10) + 0.5)))
    return(trupry1)}
  
  if(r==81) {
    trupry1 = (((X -1.249999985)^2 / 7.49999985) + 0.691666668) + ((((Z + .125)^2 / 1.25) - .0125) *((((X - 1)^2 / 2) + 0.3) - (((X -1.249999985)^2 / 7.49999985) + 0.691666668)))
    return(trupry1)}
  
  if(r==90) {
    trupry1 = (((X -1.253012048)^2 / -3.012048193) + 0.921253012) + ((((Z - 1.75)^2 / -2.5) + 1.225) *((((X -1.257575758)^2 / -7.575757576) + 0.508757576) - (((X -1.253012048)^2 / -3.012048193) + 0.921253012)))
    return(trupry1)}
  
  if(r==91) {
    trupry1 = (((X - 1.25)^2 / -2.5) + 0.925) + ((((Z + .75)^2 / 2.5) - .225) *((((X - 1.25)^2 / -5) + 0.5125) - (((X - 1.25)^2 / -2.5) + 0.925)))
    return(trupry1)}
  
  if(r==84) {
    trupry1 = (((X -1.247863248)^2 / 2.136752137) + 0.071247863) + ((((Z - 1.75)^2 / -2.5) + 1.225)*((((X -1.246268657)^2 / 3.731343284) + 0.483746269) - (((X -1.247863248)^2 / 2.136752137) + 0.071247863)))
    return(trupry1)}
  
  if(r==98) {
    trupry1 = (((X - 1)^2 / 10) + 0.8) + ((((Z + .125)^2 / 1.25) - .0125)*((((X -1.25030012)^2 / 3.00120048) + 0.07912497) - (((X - 1)^2 / 10) + 0.8)))
    return(trupry1)}
  
  if(r==99) {
    trupry1 = ((1-Z)*((((X - 1)^2 / -2) + 0.8))) + ((Z)*(((X - 1)^2 / -10 + 0.3)))
    return(trupry1)}
  
  if(r==102) {
    trupry1 = ((1-Z)*((((X + 0)^2 / 1.666666667) + 0.1))) + ((Z)*(((X + 0.250750751)^2 / 7.507507508 + 0.691624925)))
    return(trupry1)}
  
  if(r==103) {
    trupry1 = (((X + 0.249785775)^2 / 2.142245073) + 0.070874979) + ((((Z - 1.125)^2 / -1.25) + 1.0125) *((((X + 0.25)^2 / 5) + 0.5875) - (((X + 0.249785775)^2 / 2.142245073) + 0.070874979)))
    return(trupry1)}
  
  if(r==113) {
    trupry1 = ((1-Z)*((((X -5.55112E-16)^2 / -10) + 0.3))) + ((Z)*(((X + 0)^2 / -1.666666667 + 0.9)))
    return(trupry1)}
  
  if(r==108) {
    trupry1 = (((X + 0)^2 / -1.428571429) + 0.9) + ((((Z + .75)^2 / 2.5) - .225)*((((X + 0.257575758)^2 / -7.575757576) + 0.308757576) - (((X + 0)^2 / -1.428571429) + 0.9)))
    return(trupry1)}
  
  if(r==77) {
    trupry1 = 0.2+(0.1*X)+(0.2*Z)+(0.3*X*Z)
    return(trupry1)}
  
  if(r==74) {
    trupry1 = 0.2+(-0.1*X)+(0.4*Z)+(-0.3*X*Z)
    return(trupry1)}
  
  if(r==75) {
    trupry1 = 0.2+(0.2*X)+(0.1*Z)+(0.3*X*Z)
    return(trupry1)}
  
  if(r==76) {
    trupry1 = 0.6+(0.2*X)+(-0.5*Z)+(0.3*X*Z)
    return(trupry1)}
  
  if(r==78) {
    trupry1 = 0.1+(0.6*X)+(0.5*Z)+(-0.3*X*Z)
    return(trupry1)}
  
  if(r==79) {
    trupry1 = 0.1+(0.7*X)+(0.4*Z)+(-0.3*X*Z)
    return(trupry1)}
  
  if(r==93) {
    trupry1 = 0.4+(-0.1*X)+(0.5*Z)+(-0.4*X*Z)
    return(trupry1)}
  
  if(r==94) {
    trupry1 = 0.7+(-0.5*X)+(0.2*Z)+(0.4*X*Z)
    return(trupry1)}
  
  if(r==95) {
    trupry1 = 0.3+(0.6*X)+(-0.1*Z)+(-0.4*X*Z)
    return(trupry1)}
  
  if(r==92) {
    trupry1 = 0.9+(-0.3*X)+(-0.1*Z)+(-0.4*X*Z)
    return(trupry1)}
  
  if(r==105) {
    trupry1 = 0.8+(-0.6*X)+(0.1*Z)+(0.5*X*Z)
    return(trupry1)}
  
  if(r==104) {
    trupry1 = 0.9+(-0.7*X)+(-0.6*Z)+(0.5*X*Z)
    return(trupry1)}
  
  if(r==48) {
    trupry1 = (((X + 0)^2 / 3.333333333) + 0.3) + ((((Z + .75)^2 / 2.5) - .225)*((((X + 3.88578E-16)^2 / 3.333333333) + 0.6) - (((X + 0)^2 / 3.333333333) + 0.3)))
    return(trupry1)}
  
  if(r==40) {
    trupry1 = (((X - 1)^2 / 3.333333333) + 0.5) + ((((Z + .125)^2 / 1.25) - .0125)*((((X - 1)^2 / 3.333333333) + 0.2) - (((X - 1)^2 / 3.333333333) + 0.5)))
    return(trupry1)}
  
  if(r==41) {
    trupry1 = (((X - 1.25)^2 / 5) + 0.3875) + ((((Z + .125)^2 / 1.25) - .0125)*((((X - 1.25)^2 / 5) + 0.0875) - (((X - 1.25)^2 / 5) + 0.3875)))
    return(trupry1)}
  
  if(r==43) {
    trupry1 = (Z*((((X + 0)^2 / -3.333333333) + 0.4))) + ((1-Z)*(((X + 3.88578E-16)^2 / -3.333333333 + 0.8)))
    return(trupry1)}
  
  if(r==42) {
    trupry1 = (((X + 0.25)^2 / -5) + 0.9125) + ((((Z - 1.125)^2 / -1.25) + 1.0125)*((((X + 0.25)^2 / -5) + 0.5125) - (((X + 0.25)^2 / -5) + 0.9125)))
    return(trupry1)}
  
  if(r==46) {
    trupry1 = (Z*((((X + 3.88578E-16)^2 / 3.333333333) + 0.6))) + ((1-Z)*(((X + 0)^2 / 3.333333333 + 0.1)))
    return(trupry1)}
  
  if(r==49) {
    trupry1 = (((X - 1)^2 / -2.5) + 0.5) + ((((Z - 1.125)^2 / -1.25) + 1.0125)*((((X - 1)^2 / -2.5) + 0.8) - (((X - 1)^2 / -2.5) + 0.5)))
    return(trupry1)}
  
  if(r==47) {
    trupry1 = (((X -1.249625187)^2 / -3.748125937) + 0.916625037) + ((((Z - 1.75)^2 / -2.5) + 1.225)*((((X -1.249625187)^2 / -3.748125937) + 0.616625037) - (((X -1.249625187)^2 / -3.748125937) + 0.916625037)))
    return(trupry1)}
  
  if(r==44) {
    trupry1 = (((X -1.249625187)^2 / 3.748125937) + 0.083374963) + ((((Z + .75)^2 / 2.5) - .225)*((((X -1.249625187)^2 / 3.748125937) + 0.483374963) - (((X -1.249625187)^2 / 3.748125937) + 0.083374963)))
    return(trupry1)}
  
  if(r==45) {
    trupry1 = (((X + 0.25030012)^2 / -3.00120048) + 0.62087503) + ((((Z - 1.75)^2 / -2.5) + 1.225)*((((X + 0.25030012)^2 / -3.00120048) + 0.92087503) - (((X + 0.25030012)^2 / -3.00120048) + 0.62087503)))
    return(trupry1)}
  
  if(r==17) {
    trupry1 = 0.4+(-0.3*X)+(0.3*Z)+(0*X*Z)
    return(trupry1)}
  
  if(r==18) {
    trupry1 = 0.5+(-0.3*X)+(0.3*Z)+(0*X*Z)
    return(trupry1)}
  
  if(r==19) {
    trupry1 = 0.6+(0.3*X)+(-0.3*Z)+(0*X*Z)
    return(trupry1)}
  
  if(r==22) {
    trupry1 = 0.2+(0.3*X)+(0.4*Z)+(0*X*Z)
    return(trupry1)}
  
  if(r==20) {
    trupry1 = 0.6+(0.3*X)+(-0.4*Z)+(0*X*Z)
    return(trupry1)}
  
  if(r==23) {
    trupry1 = 0.1+(0.3*X)+(0.5*Z)+(0*X*Z)
    return(trupry1)}
  
  if(r==21) {
    trupry1 = 0.5+(0.4*X)+(-0.3*Z)+(0*X*Z)
    return(trupry1)}
  
  if(r==15) {
    trupry1 = 0.9+(-0.4*X)+(-0.3*Z)+(0*X*Z)
    return(trupry1)}
  
  if(r==24) {
    trupry1 = 0.1+(0.4*X)+(0.4*Z)+(0*X*Z)
    return(trupry1)}
  
  if(r==16) {
    trupry1 = 0.9+(-0.5*X)+(-0.3*Z)+(0*X*Z)
    return(trupry1)}
  
  if(r==59) {
    trupry1 = (((X+0)^2 / 3.333333) + 0.3) + ((((Z + .75)^2 / 2.5) - .225) * ((((X + .249962502)^2 / 3.749812509) + 0.4833375) - (((X + 0)^2 / 3.333333) + 0.3)))
    return(trupry1)}
  
  if(r==54) {
    trupry1 = (((X + -1)^2 / 3.33333333) + 0.5) + ((((Z + .125)^2 / 1.25) - .0125) * ((((X + -1)^2 / 5) + 0.2) - (((X + -1)^2 / 3.33333333) + 0.5)))
    return(trupry1)}
  
  if(r==55) {
    trupry1 = (((X + -1.25)^2 / 5) + 0.3875) + ((((Z + 0.125)^2 / 1.25) - 0.0125) * ((((X + -1)^2 / 5) + 0.2) - (((X + -1.25)^2 / 5) + 0.3875)))
    return(trupry1)}
  
  if(r==56) {
    trupry1 = (Z*(((X + 0)^2 / -3.33333333) + 0.4)) + ((1-Z)*(((X + 4.44089e-16)^2 / -2.5) + 0.8))
    return(trupry1)}
  
  if(r==60) {
    trupry1 = (Z*((((X + .25)^2 / 5) + 0.5875))) + ((1-Z)*(((X + 0)^2 / 2.5 + 0.1)))
    return(trupry1)}
  
  if(r==61) {
    trupry1 = (((X - 1)^2 / -2.5) + 0.5) + ((((Z - 1.125)^2 / -1.25) + 1.0125)*((((X - 1)^2 / -3.33333333) + 0.8) - (((X - 1)^2 / -2.5) + 0.5)))
    return(trupry1)}
  
  if(r==58) {
    trupry1 = (((X -1)^2 / -2) + 0.9) + ((((Z - 1.75)^2 / -2.5) + 1.225)*((((X -1.249625187)^2 / -3.748125937) + 0.616625037) - (((X -1)^2 / -2) + 0.9)))
    return(trupry1)}
  
  if(r==57) {
    trupry1 = (((X + 0.25)^2 / -2.5) + 0.725) + ((((Z - 1.75)^2 / -2.5) + 1.225)*((((X + 0.25030012)^2 / -3.00120048) + 0.92087503) - (((X + 0.25)^2 / -2.5) + 0.725)))
    return(trupry1)}
  
  if(r==51) {
    trupry1 = 0.5+(-0.3*X)+(0.3*Z)+(.1*X*Z)
    return(trupry1)}
  
  if(r==52){
    trupry1 = 0.6+(0.3*X)+(-0.4*Z)+(0.1*X*Z)
    return(trupry1)}
  
  if(r==53){
    trupry1 = 0.1+(0.3*X)+(0.5*Z)+(-0.1*X*Z)
    return(trupry1)}
  
  if(r==50){
    trupry1 = 0.9+(-0.4*X)+(-0.3*Z)+(0.1*X*Z)
    return(trupry1)}
  
  if(r==66){
    trupry1 = (((X - 1)^2 / 5) + 0.3) + ((((Z - 1.125)^2 / -1.25) + 1.0125) *((((X -1)^2 / 2.5) + .5) - (((X -1)^2 / 5) + .3)))
    return(trupry1)}
  
  if(r==69){
    trupry1 = (((X + 0)^2 / 2.5) + 0.5) + ((((Z + .125)^2 / 1.25) - .0125) *((((X + 0)^2 / 5) + 0.1) - (((X + 0)^2 / 2.5) + 0.5)))
    return(trupry1)}
  
  if(r==70){
    trupry1 = ((1-Z)*((((X + 0.2499994)^2 / 2.9999976) + 0.27916675))) + ((Z)*(((X + 0)^2 / 3.33333333 + 0.1)))
    return(trupry1)}
  
  if(r==71){
    trupry1 = (((X - 1.25)^2 / -5) + 0.9125) + ((((Z - 1.125)^2 / -1.25) + 1.0125) *((((X - 1.25030012)^2 / -3.00120048) + 0.82087503) - (((X - 1.25)^2 / -5) + 0.9125)))
    return(trupry1)}
  
  if(r==72){
    trupry1 = (((X + 0)^2 / 10) + 0.7) + ((((Z + .125)^2 / 1.25) - .0125)*((((X + 0.25)^2 / 5) + 0.1875) - (((X + 0)^2 / 10) + 0.7)))
    return(trupry1)}
  
  if(r==67){
    trupry1 = (((X + 0)^2 / -5) + 0.5) + ((((Z - 1.75)^2 / -2.5) + 1.225)*((((X + 0)^2 / -2.5) + 0.9) - (((X + 0)^2 / -5) + 0.5)))
    return(trupry1)}
  
  if(r==73){
    trupry1 = (((X -1.253012048)^2 / -3.012048193) + 0.921253012) + ((((Z - 1.75)^2 / -2.5) + 1.225) *((((X -1.25)^2 / -5) + 0.6125) - (((X -1.253012048)^2 / -3.012048193) + 0.921253012)))
    return(trupry1)}
  
  if(r==68){
    trupry1 = (((X -1.247863248)^2 / 2.136752137) + 0.071247863) + ((((Z - 1.75)^2 / -2.5) + 1.225)*((((X -1.25030012)^2 / 3.00120048) + 0.37912497) - (((X -1.247863248)^2 / 2.136752137) + 0.071247863)))
    return(trupry1)}
  
  if(r==64){
    trupry1 = 0.2+(0.1*X)+(0.3*Z)+(0.2*X*Z)
    return(trupry1)}
  
  if(r==62){
    trupry1 = 0.2+(-0.1*X)+(0.3*Z)+(-0.2*X*Z)
    return(trupry1)}
  
  if(r==65){
    trupry1 = 0.2+(0.2*X)+(0.1*Z)+(0.2*X*Z)
    return(trupry1)}
  
  if(r==63){
    trupry1 = 0.5+(0.2*X)+(-0.4*Z)+(0.2*X*Z)
    return(trupry1)}
  
  if(r==115){
    trupry1 = (((X -1)^2 / -1.66666667) + 0.8) + ((((Z + .125)^2 / 1.25) - .0125) *((((X -1.248502994)^2 / -14.97005988) + 0.90412515) - (((X -1)^2 / -1.66666667) + 0.8)))
    return(trupry1)}
  
  if(r==109){
    trupry1 = (((X -1)^2 / 1.428571429) + 0.2) + ((((Z - 1.75)^2 / -2.5) + 1.225) *((((X - 1.250750751)^2 / 7.507507508) + 0.091624925) - (((X -1)^2 / 1.428571429) + 0.2)))
    return(trupry1)}
  
  if(r==110){
    trupry1 = (((X - 1)^2 / 10) + 0.8) + ((((Z + .125)^2 / 1.25) - .0125)*((((X -1.25)^2 / 2.5) + 0.075) - (((X - 1)^2 / 10) + 0.8)))
    return(trupry1)}
  
  if(r==114){
    trupry1 = ((1-Z)*((((X - 1)^2 / -1.66666667) + 0.9))) + ((Z)*(((X - 1)^2 / -10 + 0.3)))
    return(trupry1)}
  
  if(r==106){
    trupry1 = 0.7+(-0.6*X)+(0.2*Z)+(0.5*X*Z)
    return(trupry1)}
  
  if(r==107){
    trupry1 = 0.3+(0.6*X)+(-0.1*Z)+(-0.5*X*Z)
    return(trupry1)}
  
  
  if(r==1){
    trupry1=0.1+0*X+0*Z
    return(trupry1)}
  
  if(r==116){
    trupry1=(.1+.1)-(Z*.1)+0*X
    return(trupry1)}
  
  if(r==117){
    trupry1=.1+(Z*.2)+0*X
    return(trupry1)}
  
  if(r==118){
    trupry1=.1+(Z*.3)+0*X
    return(trupry1)}
  
  if(r==13){
    trupry1=.1+(Z*.4)+0*X
    return(trupry1)}
  
  if(r==119){
    trupry1=(.1+.5)-(Z*.5)+0*X
    return(trupry1)}
  
  if(r==2){
    trupry1=0.2+0*X+0*Z
    return(trupry1)}
  
  if(r==120){
    trupry1=(.2+.1)-(Z*.1)+0*X
    return(trupry1)}
  
  if(r==121){
    trupry1=(.2+.2)-(Z*.2)+0*X
    return(trupry1)}
  
  if(r==14){
    trupry1=.2+(Z*.3)+0*X
    return(trupry1)}
  
  if(r==122){
    trupry1=.2+(Z*.4)+0*X
    return(trupry1)}
  
  if(r==123){
    trupry1=(.2+.5)-(Z*.5)+0*X
    return(trupry1)}
  
  if(r==3){
    trupry1=0.3+0*X+0*Z
    return(trupry1)}
  
  if(r==124){
    trupry1=.3+(Z*.1)+0*X
    return(trupry1)}
  
  if(r==125){
    trupry1=.3+(Z*.2)+0*X
    return(trupry1)}
  
  if(r==126){
    trupry1=.3+(Z*.3)+0*X
    return(trupry1)}
  
  if(r==10){
    trupry1=(.3+.4)-(Z*.4)+0*X
    return(trupry1)}
  
  if(r==127){
    trupry1=(.3+.5)-(Z*.5)+0*X
    return(trupry1)}
  
  if(r==4){
    trupry1=0.4+0*X+0*Z
    return(trupry1)}
  
  if(r==128){
    trupry1=(.4+.1)-(Z*.1)+0*X
    return(trupry1)}
  
  if(r==129){
    trupry1=.4+(Z*.2)+0*X
    return(trupry1)}
  
  if(r==130){
    trupry1=.4+(Z*.3)+0*X
    return(trupry1)}
  
  if(r==11){
    trupry1=(.4+.4)-(Z*.4)+0*X
    return(trupry1)}
  
  if(r==131){
    trupry1=.4+(Z*.5)+0*X
    return(trupry1)}
  
  if(r==5){
    trupry1=0.5+0*X+0*Z
    return(trupry1)}
  
  if(r==12){
    trupry1=(.5+.1)-(Z*.1)+0*X
    return(trupry1)}
  
  if(r==132){
    trupry1=(.5+.2)-(Z*.2)+0*X
    return(trupry1)}
  
  if(r==133){
    trupry1=(.5+.3)-(Z*.3)+0*X
    return(trupry1)}
  
  if(r==134){
    trupry1=.5+(Z*.4)+0*X
    return(trupry1)}
  
  if(r==6){
    trupry1=0.6+0*X+0*Z
    return(trupry1)}
  
  if(r==135){
    trupry1=.6+(Z*.1)+0*X
    return(trupry1)}
  
  if(r==136){
    trupry1=(.6+.2)-(Z*.2)+0*X
    return(trupry1)}
  
  if(r==137){
    trupry1=.6+(Z*.3)+0*X
    return(trupry1)}
  
  if(r==7){
    trupry1=0.7+0*X+0*Z
    return(trupry1)}
  
  if(r==138){
    trupry1=.7+(Z*.1)+0*X
    return(trupry1)}
  
  if(r==139){
    trupry1=(.7+.2)-(Z*.2)+0*X
    return(trupry1)}
  
  if(r==8){
    trupry1=0.8+0*X+0*Z
    return(trupry1)}
  
  if(r==140){
    trupry1=(.8+.1)-(Z*.1)+0*X
    return(trupry1)}
  
  if(r==9){
    trupry1=0.9+0*X+0*Z
    return(trupry1)}
  
  if(r==25){
    trupry1=(.1+.1)-(((((Z-1.125)^2)/-1.25)+1.0125)*.1)+0*X
    return(trupry1)}
  
  if(r==141){
    trupry1=.1+(((((Z+.125)^2)/1.25)-.0125)*.2)+0*X
    return(trupry1)}
  
  if(r==33){
    trupry1=.1+(((((Z+.125)^2)/1.25)-.0125)*.3)+0*X
    return(trupry1)}
  
  if(r==142){
    trupry1=.1+(((((Z+.75)^2)/2.5)-.225)*.4)+0*X
    return(trupry1)}
  
  if(r==143){
    trupry1=(.1+.5)-(((((Z-1.75)^2)/-2.5)+1.225)*.5)+0*X
    return(trupry1)}
  
  if(r==26){
    trupry1=(.2+.1)-(((((Z-1.75)^2)/-2.5)+1.225)*.1)+0*X
    return(trupry1)}
  
  if(r==27){
    trupry1=(.2+.2)-(((((Z-1.75)^2)/-2.5)+1.225)*.2)+0*X
    return(trupry1)}
  
  if(r==144){
    trupry1=.2+(((((Z+.75)^2)/2.5)-.225)*.3)+0*X
    return(trupry1)}
  
  if(r==145){
    trupry1=.2+(((((Z+.125)^2)/1.25)-.0125)*.4)+0*X
    return(trupry1)}
  
  if(r==28){
    trupry1=(.2+.5)-(((((Z-1.75)^2)/-2.5)+1.225)*.5)+0*X
    return(trupry1)}
  
  if(r==146){
    trupry1=.3+(((((Z+.125)^2)/1.25)-.0125)*.1)+0*X
    return(trupry1)}
  
  if(r==147){
    trupry1=.3+(((((Z+.125)^2)/1.25)-.0125)*.2)+0*X
    return(trupry1)}
  
  if(r==34){
    trupry1=.3+(((((Z+.125)^2)/1.25)-.0125)*.3)+0*X
    return(trupry1)}
  
  if(r==148){
    trupry1=(.3+.4)-(((((Z-1.125)^2)/-1.25)+1.0125)*.4)+0*X
    return(trupry1)}
  
  if(r==29){
    trupry1=(.3+.5)-(((((Z+.75)^2)/2.5)-.225)*.5)+0*X
    return(trupry1)}
  
  if(r==30){
    trupry1=(.4+.1)-(((((Z-1.125)^2)/-1.25)+1.0125)*.1)+0*X
    return(trupry1)}
  
  if(r==35){
    trupry1=.4+(((((Z+.75)^2)/2.5)-.225)*.2)+0*X
    return(trupry1)}
  
  if(r==149){
    trupry1=.4+(((((Z-1.125)^2)/-1.25)+1.0125)*.3)+0*X
    return(trupry1)}
  
  if(r==150){
    trupry1=(.4+.4)-(((((Z+.75)^2)/2.5)-.225)*.4)+0*X
    return(trupry1)}
  
  if(r==36){
    trupry1=.4+(((((Z-1.75)^2)/-2.5)+1.225)*.5)+0*X
    return(trupry1)}
  
  if(r==151){
    trupry1=(.5+.1)-(((((Z+.75)^2)/2.5)-.225)*.1)+0*X
    return(trupry1)}
  
  if(r==31){
    trupry1=(.5+.2)-(((((Z+.75)^2)/2.5)-.225)*.2)+0*X
    return(trupry1)}
  
  if(r==152){
    trupry1=(.5+.3)-(((((Z+.75)^2)/2.5)-.225)*.3)+0*X
    return(trupry1)}
  
  if(r==37){
    trupry1=.5+(((((Z-1.125)^2)/-1.25)+1.0125)*.4)+0*X
    return(trupry1)}
  
  if(r==153){
    trupry1=.6+(((((Z-1.125)^2)/-1.25)+1.0125)*.1)+0*X
    return(trupry1)}
  
  if(r==154){
    trupry1=(.6+.2)-(((((Z+.125)^2)/1.25)-.0125)*.2)+0*X
    return(trupry1)}
  
  if(r==38){
    trupry1=.6+(((((Z-1.125)^2)/-1.25)+1.0125)*.3)+0*X
    return(trupry1)}
  
  if(r==39){
    trupry1=.7+(((((Z-1.75)^2)/-2.5)+1.225)*.1)+0*X
    return(trupry1)}
  
  if(r==32){
    trupry1=(.7+.2)-(((((Z+.75)^2)/2.5)-.225)*.2)+0*X
    return(trupry1)}
  
  if(r==155){
    trupry1=(.8+.1)-(((((Z+.125)^2)/1.25)-.0125)*.1)+0*X
    return(trupry1)}
  
}

##########
# Draw Continuous IV Datasets
##########

# Prepare filenames
filenames<-paste("sim0", 1:9, ".dta", sep = "")
filenames<-c(filenames,paste("sim", 10:155, ".dta", sep = ""))

ychange<-function(k){ifelse(k[1]<k[2],1,0)}

for(r in 1:155){

	print(c(r, " Continuous"), quote=F, sep="")
	X<-runif(100000, min=0, max=1)
	Z<-runif(100000, min=0, max=1)
	pry1 = trupry(r,X,Z)
	temp<-runif(100000, min=0, max=1)
	y<-apply(cbind(temp, pry1), 1, ychange)
	x1<-X
	x2<-Z
	x1x2<-x2*x1
	x1sq<-x1*x1
	x2sq<-x2*x2
	x1sqx2<-x1sq*x2
	x2sqx1<-x2sq*x1
	x1sqx2sq<-x1sq*x2sq

	dat<-cbind(x1,x2,x1x2,x1sq,x2sq,x1sqx2,x2sqx1,x1sqx2sq,pry1,y)
	dat<-as.data.frame(dat)
	dimnames(dat)[[2]]<-c("x1","x2","x1x2","x1sq","x2sq","x1sqx2","x2sqx1","x1sqx2sq","pry1","y")
	library(foreign)
	write.dta(dat, filenames[r],version=7)
	
	rm(X,Z,x1,x2,pry1,temp,y,x1x2,x1sq,x2sq,x1sqx2,x2sqx1,x1sqx2sq,dat)

}

##########
# Dichotomize Continuous IV Datasets
##########

# Prepare filenames
dfilenames<-paste("sim0", 1:9, "D.dta", sep = "")
dfilenames<-c(dfilenames,paste("sim", 10:155, "D.dta", sep = ""))

for(r in 1:155){

	print(c(r, " Dichotomous"), quote=F, sep="")

	read.dta(filenames[r])->simdatmain
	
	simdatmain$x2<-round(simdatmain$x2,digits = 0)
	simdatmain$x1x2<-simdatmain$x2*simdatmain$x1
	simdatmain$x2sq<-simdatmain$x2*simdatmain$x2
	simdatmain$x1sqx2<-simdatmain$x1sq*simdatmain$x2
	simdatmain$x2sqx1<-simdatmain$x2sq*simdatmain$x1
	simdatmain$x1sqx2sq<-simdatmain$x1sq*simdatmain$x2sq
	simdatmain$pry1<-trupry(r,simdatmain$x1,simdatmain$x2)
	temp<-runif(length(simdatmain$x2))
	simdatmain$y<-apply(cbind(temp,simdatmain$pry1), 1, ychange)
	
	write.dta(simdatmain, dfilenames[r], version=7)
	rm(simdatmain)

}